#include "Image.H"
#include "IO.H"
#include "Util.H"

#ifndef graymophologyH
#define graymophologyH




void discStructureElement(pip::Image<unsigned char> &imOut, int rad);

template <class Tin>
void grayDilate(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
	//use a function pointer
	bool (*fp)(Tin,Tin)= Dilate;
	doWork(imIn,imOut,se,(*fp)); // someone else is going to do the work
}

template <class Tin>
void grayErode(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
	bool (*fp)(Tin,Tin)= Erode; // just use another function
	doWork(imIn,imOut,se,(*fp));
}

template <class Tin>
void grayOpening(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
	pip::Image<Tin> tmp;
	grayErode( imIn, tmp, se);  // its just a matter of definition ... 
	grayDilate( tmp, imOut, se);
}

template <class Tin>
void grayClosing(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
	pip::Image<Tin> tmp;
	grayDilate( imIn, tmp, se); // the other way around...
	grayErode( tmp, imOut, se);
}

// some functions to decide  wether '>' or '<' !
template <class Tin>
bool Dilate(Tin image,Tin value){
	return (image > value);
}	

template <class Tin>
bool Erode(Tin image,Tin value){
	return (image < value);
}	

//unused functions, since the gray-variants do the binary tasks as well
template <class Tin>
void Opening(pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
}

template <class Tin>
void Closing(pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se){
}

// here the work is done
template <class Tin>
void doWork(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut, const pip::Image<unsigned char> &se,bool (*fp)(Tin,Tin)){
	imOut.resize(imIn.dimension());  //prepare the output to be filled
	//init...
	Tin value(0), init_value(0);      
	int posx(0),posy(0),kern_wi(se.dim(0)/2);
	float min,max,avr;
	getMinMaxAvrImage(imIn,min,max,avr);

	//determine init_value corresponding to function used - so we can init the search every time
	if((*fp)(max,min))
		init_value=min;
	else
		init_value=max;
	//over all pixel
	for(int x=0;x<imIn.dim(0);++x){
		for(int y=0;y<imIn.dim(1);++y){
			value=init_value;                             // init the search value ...!
			//apply kernel
			for(int     ix=-kern_wi ; ix <= kern_wi ; ++ix){
				for(int iy=-kern_wi ; iy <= kern_wi ; ++iy){
					if( se(ix+kern_wi,iy+kern_wi)==1){ // search within '1' pixel in kernel
						posx = ix+x;
						posy = iy+y;
						if(	   posx > 0 
							&& posx < imIn.dim(0)
							&& posy < imIn.dim(1)
							&& posy > 0           // and make sure we are inside the image ... I assumed that borders might be illegal, as they may change values...
						   )
						   if((*fp)(imIn(posx,posy),value))   //search for 'extreme' value in 'kernel-shape'
							   value=imIn(posx,posy);		  //and remember it
					}
				imOut(x,y)=value;								// the best one belongs into the output
				}
			}
		}
	}
	// getMinMaxAvrImage(imOut,min,max,avr);  //if someone want to do a COUT or something
}




// this function is out of use ... it should have applied the operators for all levels of a 3d image... 
/*template <class Tin>
void discWise(pip::Image<Tin> &imIn,
			  const pip::Image<unsigned char> &se,
			  void (*fctPT)(pip::Image<Tin> &, const pip::Image<unsigned char> &,int)){
				  for(int z=0; z < imIn.dim(2); ++z){	
					(*fctPT)(imIn,se,z);
				  }
}*/


//some old attempts to use a UMBRA... however- my virtual machine had problems to handle 150mb images... (which is an 3d image of the fundus -already)

/*

template <class Tin>
void constructUmbra(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut){
	float min,max,avr;
	int zscale,offset,first,second;
	first= imIn.dim(0);
	second=imIn.dim(1);
	getMinMaxAvrImage(imIn,min,max,avr);
	zscale = abs(int(max)-int(min)); //how large should the z-axis be?
	offset = int(min); //need this to fit the values of inIm into new z-axis

	std::cout << first <<"  " << second <<"   " << zscale;
	pip::Image<Tin> tmp(Dimension(first,second,zscale)); // resize/create and fill the image
		
	for(int x=0; x<imIn.dim(0);++x){
		for(int y=0; y<imIn.dim(1);++y){
			for(int z=0; z<imIn.dim(2);++z){
				if(z<=(int(imIn(x,y))+offset))  //fill all beneath the value (shifted origin)
					//tmp(x,y,z)=1; 
					std::cout;
				else
					z=imIn.dim(2); //end this column
			}
		}
	}
	//imOut=tmp;
}

template <class Tin>
void collapseUmbra(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut){
	float min,max,avr;
	int zscale,offset;
	getMinMaxAvrImage(imIn,min,max,avr);
	offset = int(min); //need this to fit the values of inIm into new z-axis
	pip::Image<Tin> tmp(imIn);

	for(int x=0; x<imOut.dim(0);++x){
		for(int y=0; y<imOut.dim(1);++y){
			for(int z=0; z<imOut.dim(2);++z){
				if( 1==imOut(x,y,z)  
					&& (z==(imOut.dim(2)-1)
					    || 0==imOut(x,y,z)))
					tmp(x,y)=z-offset;
			}
		}
	}
	imOut.resize(Dimension(imIn.dim(0),imIn.dim(1))); // resize/create and fill the image
	imOut=tmp;
}

template <class Tin>
void constructTop(const pip::Image<Tin> &imIn, pip::Image<Tin> &imOut){
	float min,max,avr;
	int zscale,offset;
	getMinMaxAvrImage(imIn,min,max,avr);
	zscale = abs(int(max)-int(min)); //how large should the z-axis be?
	offset = int(min);
    

	imOut.resize(imIn.dim(0),imIn.dim(1),zscale); // resize/create and fill the image
	imOut.fill(0);

	if(!imIn.dim(2)>0){ //raw gray-scale image
		for(int x=0; x<imOut.dim(0);++x){
			for(int y=0; y<imOut.dim(1);++y){
				imOut(x,y,(int(imIn(x,y))+offset))=1; 
			}
		}
	}else{  //its an umbra
		for(int x=0; x<imOut.dim(0);++x){
			for(int y=0; y<imOut.dim(1);++y){
				for(int z=0; z<imOut.dim(2);++z){
					if(imIn(x,y,z)==0)		// if the umbra (always 1) ends ...
						imOut(x,y,z-1)=1;  //we go one step back:  top surface of umbra
				}
			}
		}
	}
}*/

	
	/*
	 OLD DILATE

	pip::Image<Tin> imIn(imOut);

	int posx(0),posy(0);
	// go through the whole slice 
	for(int _x=0;_x<imIn.dim(0);++_x){
		for(int _y=0;_y<imIn.dim(1);++_y){
			if(imIn(_x,_y,z)==1){  // if there is a pixel in the 'original' ...
				//over the 'SE' elements
				for(int ix=-(se.dim(0)/2);ix<(se.dim(0)/2);++ix){
					for(int iy=-(se.dim(1)/2);iy<(se.dim(1)/2);++iy){
						posx = ix+_x;
						posy = iy+_y;
						if(    se(ix+se.dim(0)/2,iy+se.dim(1)/2)==1
							&& posx > 0 
							&& posx < imIn.dim(0)
							&& posy < imIn.dim(1)
							&& posy > 0         
						   )
							imOut(posx,posy)=0;
					}
				}
			}
		}
	}*/



#endif

